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STOCHASTIC  SUBSURFACE  FLOW  AND  TRANSPORT  IN  FRACTAL 
CONDUCTIVITY  FIELDS 

ALBERTO  S.  NDUMU,  PAUL  S.  ADDISON 
Civil  Engineering  Group,  School  of  the  Built  Environment,  Napier  University, 

10  Colint  on  Road,  Edinburgh,  EH  10  5DT,  Scotland,  UK. 


Monte  Carlo  simulations  of  subsurface  flow  and  contaminant  transport  of  a non-reactive  solute 
plume  by  steady-state  flow  with  a uniform  velocity  were  performed  in  a two-dimensional  synthetic 
heterogeneous  porous  media  whose  hydraulic  conductivity  is  non-stationary  and  described  by  multi- 
scale fractional  Brownian  motion.  Analysis  of  the  flow  and  transport  results  indicates  that  the 
longitudinal  velocity  variance  is  nearly  constant  in  the  longitudinal  direction  while  in  the  transverse 
direction  it  assumes  a parabolic  shape.  The  velocity  variance  is  maximum  at  the  impervious 
boundaries  and  decreases  in  transverse  direction  with  distance  from  the  boundaries  reaching  the 
minimum  value  at  the  domain  centre.  We  observe  that  the  particle  displacement  covariance  is 
anomalous  or  non-Fickian  at  all  times  t in  the  dispersion  process  irrespective  of  the  Hurst  exponent 
H and  grows  temporally  faster  than  linearly. 


1 Introduction 

Natural  variability  of  subsurface  (i.e.,  hydrologic  and  geologic  systems)  hydraulic 
properties  is  the  main  factor  controlling  the  flow  and  spreading  of  contaminants  in 
subsurface  porous  media.  The  primary  physical  property  which  exhibits  large-scale 
natural  spatial  variability  is  the  hydraulic  conductivity.  One  important  source  of  error  in 
numerical  models  to  simulate  flow  and  transport  in  the  subsurface  stems  from  lack  of 
knowledge  concerning  the  spatial  variability  of  hydraulic  conductivity.  In  practice,  we 
typically  have  sparse  data  which  is  inadequate  to  completely  describe  the  spatial 
variability.  This  is  a consequence  of  the  high  data  collection  cost  required  to  fully 
characterise  a given  subsurface  porous  medium.  Hence,  the  variation  of  subsurface 
hydraulic  conductivity  cannot  generally  be  described  in  all  its  detail  and  is  therefore 
uncertain1. 

One  response  to  the  problem  of  the  uncertainty  of  subsurface  hydraulic  conductivity  has 
been  an  interest  in  stochastic  methods,  which  provide  a formal  framework  for  the 
treatment  of  uncertainty.  Spatial  variability  of  subsurface  properties  is  such  that  their 
unique,  deterministic,  description  is  not  feasible2,  and  this  is  formally  recognised  in  the 
stochastic  approach,  in  which  subsurface  hydraulic  conductivity  is  commonly  regarded  as 
a random  scalar  field  (RSF). 

Stochastic  theories  involve  the  description  of  the  local  porous  medium  structure  using  a 
statistical  model  that  requires  a small  number  of  parameters  to  be  identified  from  field 
measurements.  The  detailed  spatial  distribution  of  logconductivity  is  conventionally 
reduced  to  a few  statistical  parameters,  for  example,  covariance  function,  mean,  variance, 
and  correlation  length.  It  thus  follows  that  the  outputs  from  stochastic  models  are 
probabilistic,  characterised,  for  example,  by  the  statistical  moments  or  the  full  probability 
density  function  of  the  variable  of  interest. 

In  this  study,  we  perform  Monte  Carlo  experiments  by  generating  a large  number  of 
independent  random  conductivity  realisations  with  a fractal  semivariogram  function, 
solving  the  stochastic  boundary  value  problem  by  repetitively  solving  a set  of 
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deterministic  flow  problems,  each  of  which  is  an  equally  probable  representation  of  the 
response  of  the  real  heterogeneous  medium.  The  flow  problem  is  solved  by  employing  a 
block-centred  finite  difference  scheme  (e.g.,  Smith  and  Freeze3)  to  obtain  subsurface 
velocity  fields.  A random  walk  particle  tracking  algorithm  is  used  to  solve  the  transport 
equation  (e.g.,4'7)  using  realised  velocity  fields.  Realised  solutions  are  averaged  to  get 
fluctuating  log  conductivity  semivariograms,  velocity  variances,  particle  displacement 
covariance  and  mean  solute  concentrations. 


2 Random  Field  Model 

In  the  Monte  Carlo  simulation  of  subsurface  flow  and  transport,  the  first  step  is  the  choice 
of  the  statistical  model  that  represents  the  medium  heterogeneity,  mainly  the  log  hydraulic 
conductivity,  In  Al(x) . Most  stochastic  theoretical  and  numerical  modelling  approaches 

(see  e.g.7'9)  assume  that  the  log  hydraulic  conductivity  f(x)  = In  K(x)  where  K(x)  is  the 

local  hydraulic  conductivity  of  the  subsurface  porous  medium  and  x is  the  vector  of 
spatial  co-ordinates,  as  stationary  random  scalar  field,  normally  distributed  and 
characterised  by  a constant  mean  (y)  and  isotropic  exponential  covariance  functions  e.g., 

Cy(r ) = a\e~rlx  (with  r = |r| , where  r is  the  two-point  isotropic  lag,  a\  is  the  variance 

of  the  distribution,  and  A the  finite  integral  scale).  The  justification  for  using  this  model 
is  based  on  geostatistical  data  obtained  from  several  sites  (e.g.,10).  In  addition,  the 
macrodispersivity  measured  from  tracer  tests  agree  approximately  well  with  the  stochastic 
theory  results  based  on  the  exponential  model.  Zhan  and  Wheatcraft"  argue  that,  field 
hydraulic  conductivity  measurements  are  limited,  have  large  uncertainty  and  have  been 
carried  out  for  relatively  small  spatial  scales  (at  most  a kilometre).  If  the  measurements 
are  available  at  much  larger  spatial  scales,  conductivity  values  may  remain  correlated  at 
larger  scales  and  may  give  rise  to  non-stationary  fractal  or  self-similar  distributions12'14 
with  infinite  correlation  length.  Neuman15  proposes  a spatial  scaling  assumption  with  a 
semivariogram  function  given  in  (1)  to  describe  such  a distribution  as 

ry(r)  = + r)  - y(x)]2  j = ro(r/rfH , (1) 


where  A is  a reference  scale,  y0  is  a constant  variance  parameter,  and  the  scaling 
exponent  H is  the  Hurst  exponent16.  Equation  (1)  has  been  demonstrated  to  be  valid  for 
self-affine  stochastic  processes  over  the  broad  range  0<H<\24.  The  exponent  H is 
associated  with  the  fractal  dimension  Ds  = E + 1 - H , where  E is  the  space 
dimensionality.  For  0 < H < 1 the  variance  and  integral  scale  of  the  field  Y grow  infinitely, 
thus  describing  a medium  with  spatially  evolving  heterogeneity.  The  semivariogram, 
given  by  (1),  shows  a continuous  growth  with  distance.  The  presence  of  more  than  one 
reference  scale  of  variability  has  an  influence  on  fluid  flow  and  contaminant  migration6. 
Spatially  evolving  scale  formations  are  characterised  by  a variance  which  keeps  growing 
with  the  domain  size.  For  a bounded  domain,  the  order  of  magnitude  of  the  variance  is 
given  from  (1)  by  the  following: 
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cr2y=roR2H,  (2) 

where  R is  the  characteristic  dimension  of  the  domain  (upper  fractal  cut-off  scale). 
Hence,  the  log  conductivity  integral  scale  X shows  the  following  scaling  property: 

X oc  aR  (3) 

where  a is  a constant  that  depends  on  H . For  independent  realisations  of  the  log 
conductivity  field,  the  variance  of  the  random  processes  varies  because  of  lack  of 
stationarity.  However,  averaging  over  many  realisations,  the  order  of  magnitude  of  the 
variance  is  controlled  by  the  above  relationship.  The  spatial  variability  of  the  field  is 
controlled  by  H while,  for  a given  H,  the  contrast  existing  between  the  actual  values  of 
log  conductivity  Y is  controlled  by  yQ . Increasing  y0 , keeping  H constant,  amplifies  the 
contrast  between  block  conductivites.  To  perform  hydraulic  conductivity  simulations  for 
different  H values,  the  same  order  of  magnitude  of  the  fluctuations  needs  to  be 
maintained.  This  was  needed  in  order  to  obtain  realistic  fluctuations  for  the  case  H > 1/2  . 
From  a practical  point  of  view  we  impose  in  all  cases  the  following  condition: 

y0R2H  = (7y,  (4) 


where  crj  is  the  variance  of  the  field  at  the  maximum  size  of  the  modelled  domain  (upper 
cut-off  scale).  This  means  that  once  we  fix  the  value  y0  for  the  case  H<  1/2 , say  [y^  for 
H=HX,  the  coefficient  [y0]2  for  the  case  H=H2>  1/2  becomes: 


1 Ni 

°l2  r2H2 


R2h' 


(5) 


Field  measurements  of  hydraulic  conductivity  have  indicated  an  approximately  log-normal 
distribution10.  The  same  distribution  has  therefore  been  adopted  for  the  simulations 
generated  in  this  study.  Essentially,  a normally  distributed  log  conductivity  field  Y(x,y)  of 

stationary  increments  and  isotropic  semivariogram  given  in  (1)  is  generated  over  a two- 
dimensional  domain,  the  hydraulic  conductivity  field  is  obtained  by  the  transformation 

K(x,y)  = exp[y(x,y)] , (6) 

where  K[x,y)  is  the  conductivity  assigned  to  the  point  x,y  of  the  domain.  The  turning 

band  method19'20  was  used  to  simulate  the  hydraulic  conductivity  fields.  Using  this 
approach  we  do  not  require  to  filter  the  wavelengths  larger  than  the  field  dimension 
(fractal  cut-off  scale),  from  the  Y spectrum.  The  larger  scale  of  variability  is  limited  by 
the  field  dimension,  and  the  cut-off  is  introduced  by  the  fact  that  the  generation  is 
performed  in  a finite  domain  similar  to  the  work  of  Beilin  et  al.,6.  Since  the  log 
conductivity  field  Y is  non-stationary,  the  mean  value  of  Y is  maintained  constant  in  each 
realisation  of  the  field  Y,  by  conditioning  the  field  to  a given  constant  value21.  This 
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translates  to  conditioning  the  mean  velocity  field  to  a given  constant  value.  Figure  (la) 
and  (lb)  shows  two  realised  fractal  hydraulic  conductivity  fields  generated  using  the 
turning  bands  method.  The  simulated  fields  pertain  to  H = 0.3,  and  H = 0.8  respectively. 
These  fields  show  dramatically  different  behaviour.  When  H>  1/2,  Y is  positively 
correlated  and  the  fields  shows  smooth  variations  and  large-scale  persistence22  of  the 
positive  and  negative  values,  Fig.  lb.  Conductivity  values  tend  to  cluster  above  or  below 
the  mean  for  quite  some  distance  before  they  change  to  the  other  side  of  the  mean.  When 
H < 1/2 , Y is  negatively  correlated  and  leads  to  less  persistence  within  the  conductivity 
values,  and  hence  Y varies  erratically,  Fig. la. 


(a)  H = 0.3  (b)  H = 0.8 

Figure  1:  Fractal  Hydraulic  Conductivity  Fields  Simulated  by  Turning  Bands  Method. 

To  address  model  accuracy,  we  consider  the  following  statistic: 


2 

(7) 

where  N denotes  the  number  of  pairs  of  log  hydraulic  conductivity  Y with  separation 
equal  to  r . yy(r)  is  the  semivariogram  obtained  by  spatial  averaging  all  the  equi-distant 

pairs  in  a single  realisation  over  a single  replicate  and  assuming  stationarity.  The 
reconstructed  expected  value  of  the  log  conductivity  semivariograms  for  H = 0.3, 0.8,  are 
shown  in  Figure  (2a)  and  (2b)  respectively.  The  expected  values  are  computed  in  a Monte 
Carlo  sense  averaging  over  a number  of  MC  independent  realisations  as 

Y^c^h{rm)=~frf\)  (8) 


according  to  the  following  convergence  criteria6 


n 


<s. 


(9) 


159 


where  rk  is  the  distance  of  the  two  point  lag  along  the  North-South  (NS)  and  East-West 
(EW)  directions  of  the  modelled  domain  and  M is  the  number  of  points  used  to  discretise 
the  semivariogram.  The  simulated  and  theoretical  semivariograms  (Eq.  (1))  along  NS 
(squares)  and  EW  (triangles)  directions  for  different  H values  agree  well  in  both  cases  in 
Fig.  2 although  the  accuracy  deteriorates  with  lag.  The  plotted  semivariograms  shows  that 
the  variance  grows  infinitely  with  lag  distance  for  a medium  described  by  a semivariogram 
given  in  Eq.  1 . 
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(a)  H = 0.3  (b)  H = 0.8 


Figure  2:  Simulated  and  theoretical  yy(r)  = a(r/XfH  semivariograms  alongNS  and  EW  for  fields 

in  Figure  (1). 


3 Stochastic  Finite  Difference  Flow  Model 

Steady-state  confined  flow  of  groundwater  in  a two-dimensional,  horizontal,  saturated, 
imcompressible  porous  medium  with  physical  heterogeneity  represented  by  simulated 
spatially  varying  fractal  conductivity  fields  is  considered  in  this  study.  Flow  simulations 
were  performed  in  a Monte  Carlo  manner  for  fractal  fields  of  various  H values.  The 
simulation  domain  is  a square  domain  where  constant  hydraulic  heads  are  assigned  at  the 
left  and  right  boundary  and  the  no-flow  condition  on  the  top  and  bottom  to  create  a 
uniform  mean  velocity  v . The  numerical  simulations  are  performed  for  fixed  domain  of 
252  x 252 . Beilin  el  al.3 * * *  7 suggest  that  flow  and  transport  solutions  are  not  affected  by 
refinements  of  the  grid  size  involving  more  than  four  points  per  integral  scale  2 . Here,  we 
choose  the  following  grid  sizes  in  the  longitudinal  and  transverse  directions, 
Ax  = Ay  = 2/4  . 

Steady-state  subsurface  flow  in  the  2-d  domain  is  described  by 

V.q(x)  = V.(-K(x)V^(x))  = 0 (10) 

where  q(x)  is  the  Darcy  flux  relative  to  the  solid  matrix  and  <z>(x)  is  the  hydraulic  head. 

The  flow  equation  is  discretised  by  employing  a mesh  centered-finite  difference  scheme. 
The  hydraulic  heads  at  the  nodal  points  are  solved  with  the  boundary  conditions  given 
above  by  a Gauss-Seidel  method  with  successive  over  relaxation  (SOR).  In  order  to 
compute  nodal  heads,  the  interblock  hydraulic  conductivities  are  computed  by  geometric 
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averaging  of  adjacent  conductivity  values.  This  form  of  averaging  ensures  continuity  of 
the  head  field  and  conservation  of  mass  flux  across  block  boundaries17.  Darcy’s  law  is 
applied  to  obtain  the  flux  components  q(x)  and  nodal  velocities  are  computed  by  finite 
differencing  with  a constant  porosity  value  of  tj  = 0.3 . 

In  order  to  check  the  accuracy  and  efficiency  of  the  finite  difference  approximation  of 
the  flow  equation,  we  calculated  the  local  mass  balance  error  which  can  be  defined18  as 

e = {Qou,-Qin)/Q„x  100%  (11) 

where  Qlul  and  Qm  are  the  total  volume  of  flow  out  of  the  right  boundary  and  into  the  left 
boundary,  respectively.  The  mass  balance  error  decreases  when  a smaller  value  of 
convergence  criterion,  // , in  the  iterative  scheme  is  used  (the  head  changes  in  all  nodes 
between  two  iterations  is  less  than  or  equal  to  /j  , iteration  stops).  We  found  that  for  each 
Monte  Carlo  simulation  e is  always  less  than  4%  for  each  Hurst  exponent  H . These 
result,  show  that  the  finite  difference  solution  is  reasonably  accurate,  and  accordingly,  the 
obtained  velocity  field  is  sufficiently  accurate  to  be  used  in  the  transport  model.  The  mean 
velocity  computed  for  Hurst  exponent  and  with  reference  to  1500  Monte  Carlo 
realisations  was  equal  to  unity  leading  to  the  conclusion  that  the  mean  velocity  can  indeed 
be  assumed  constant  and  unit  in  each  realisation. 


4 Stochastic  Solute  Transport  Model 

The  random  walk  particle  tracking  method  is  commonly  used  in  the  field  of  statistical 
physics  to  model  processes  involving  diffusion.  This  approach  has  also  been  used 
successfully  to  simulate  reactive  and  non-reactive  transport  in  the  subsurface4'7.  The  basic 
idea  is  to  approximate  the  spatial  distribution  of  a transport  quantity  by  a set  of  moving 
particles.  The  spatial  location  of  particles  are  updated  at  each  time  step  according  to  the 
following  equation 


X(/  + At)  = X(t)  + [v(X,t)  + V-d(v(X,t))]At  + [2d(v(X,/))A/]1/2  Z (12) 


where  x(t  + At)  is  the  updated  position  of  the  particle  that  was  at  X(z)  in  the  previous 
step,  V(x,t)  is  the  velocity  vector  at  the  old  position  at  timet,  d is  the  local  scale 

dispersion  tensor,  At  is  the  time  step,  and  Z is  a vector  of  normally  distributed  random 
numbers  of  zero  mean  and  unit  variance.  The  second  term  on  the  right-hand  side  moves 
the  particle  advectively  on  the  basis  of  the  local  velocity  field  at  each  point.  The  third  term 
is  important  when  stagnation  regions  or  sharp  fronts  exist  within  the  field23.  The  last  term 
accounts  for  the  local  scale  dispersion.  The  particle  velocity  needed  in  (12)  is  obtained  by 
using  a bilinear  interpolation  utilising  the  four  velocity  values  surrounding  the  particle 
position.  Eq.  12  is  used  to  track  2000  non-interacting  particles  initially  distributed  along  a 
strip  of  length  6X  normal  to  the  mean  flow  direction  of  the  domain.  A constant 
dimensionless  time  step  At  = 0.05  was  chosen  such  that  the  fraction  of  the  cell’s  length 
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traversed  by  a particle  in  a single  time  step  is  much  smaller  than  one.  The  dimensionless 
time  is  given  by  t = vt/X  . The  particle  tracking  experiments  are  performed  in  the  inner 
core  of  the  flow  domain  that  is  not  affected  by  the  boundaries.  Particle  tracking  is 
terminated  during  each  Monte  Carlo  run  before  the  contaminant  plume  exits  the  inner  core 
of  the  modelled  domain. 


5 Results  from  Flow  and  Transport  Stochastic  Analyses 

Numerical  generation  of  synthetic  hydraulic  conductivity  fields,  flow  computations  and 
transport  simulations  were  carried  out  in  a Monte  Carlo  manner  for  each  Hurst  exponent 
H . The  simulation  was  terminated  when  the  convergence  criteria  given  in  (9)  were 
satisfied.  For  each  Hurst  exponent  at  least  1500  Monte  Carlo  runs  were  required  to  attain 
convergence.  Numerical  analyses  were  then  carried  out  to  calculate  mean  statistical 
quantities  of  interest  from  the  flow  and  transport  simulations. 

Beilin  et  al?  show  that  the  presence  of  boundaries  influences  the  hydaulic  head  and 
velocity  variabilities  in  two-dimensional  domains.  In  Fig.  3,  we  plot  results  for  the  head 
variances  cr2H  for  different  values  of  H averaged  over  1500  Monte  Carlo  realisations. 
The  head  variability  assumes  a parabolic  shape  for  each  value  of  H . The  parabolic  shape 
of  the  curves  reflects  the  statistical  heterogeneity  of  the  standard  deviation  in  hydraulic 
head  due  to  the  constant  head  values  specified  on  the  boundaries  of  the  flow  domain.  The 
head  variance  is  zero  at  the  fixed  boundaries  as  expected  and  increases  at  the  centre  of  the 
domain  where  it  is  maximum.  The  variance  of  hydraulic  head  crj, , increases  with  H . 


Figure  3:  Head  Variance  <t2h  along  the  Longitudinal  Axis  for  different  values  of  H . 


Figs.  4a  and  4b  show  the  average  Monte  Carlo  longitudinal  velocity  variance  <j\  along 
the  center  line  in  the  longitudinal  direction  (dashed  curve),  and  also  the  center  line  along 
the  transverse  direction  (solid  curve),  for  H = 0.4  and  0.8  respectively. 
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(a)  (b) 

Figure  4:  Velocity  Variance  ay  Along  the  Centreline  in  the  Longitudinal  (Dashed)  and 
Transverse  (Solid)  Direction. 

While  in  the  longitudinal  direction  ay  is  nearly  constant,  in  the  transverse  direction  it 

assumes  a parabolic  shape.  The  longitudinal  velocity  variance  is  maximum  at  the 
impervious  boundaries  and  decreases  in  transverse  direction  with  distance  from  the 
boundaries  reaching  the  minimum  value  at  the  domain  centre. 

Travel  time  statistics  are  important7  because  they  are  robust  in  characterising  the 
dispersion  process  blending  all  sources  of  uncertainty  into  a unique  curve.  Fig.  5 
illustrates  travel  time  statistical  distributions  (breakthrough  curves)  at  three  distinct 
absorbing  barriers  placed  normal  to  the  mean  flow  respectively  at  distances  52 , 102. , and 
152  computed  by  counting  the  number  of  particles  the  cross  each  barrier  for  each 
transport  simulation  time.  The  figure  is  the  average  of  1500  Monte  Carlo  simulations  for 
the  case  H-  0.7 . This  result  shows  that  as  the  plume  travels  further  from  the  contaminant 
source,  the  increased  spreading  of  the  plume  results  in  both  a corresponding  attenuation 
and  reduction  in  peak  value  of  the  breakthrough  curves. 
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Figure  5:Travel  Time  Statistics  for  H = 0.7  . 

At  the  end  of  each  Monte  Carlo  run,  the  moments  of  the  dispersing  plume  are  computed  as 
follows:  (a)  the  trajectory  of  the  plume  centroid: 

Rr{')=-±-Ntxr-p{t), 

np  p=\ 


(13) 
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where  X",,p(t)  (i  = 1,2 ) are  the  coordinates  of  the  p -th  particle  for  the  realisation  m and 
NP  is  the  total  number  of  particles;  (b)  the  second-order  central  plume  moment: 


The  moments  (13)  and  (14)  represent  the  overall  plume  behaviour  in  each  one  of  the 
(equally  probable)  log  conductivity  fields.  Average  trajectories  (R,(i)) , and  the  average  of 

the  second  spatial  moments  are  obtained  by  taking  expectations  over  MC 

independent  Monte  Carlo  realisations: 


02468  10  02468  10 


Figure  6:  Particle  Displacement  Covariances  is 1.(0)  -S|](o)  against  t =vtjX 


for  different  values  of  H . 


The  expected  values  (15)  and  (16)  represent  the  most  probable  estimate  of  the  actual 
plume  behaviour.  Fig.  6 represents  the  average  of  1500  simulated  longitudinal  second 
spatial  moments  (511(t))-5'11(o)as  a function  of  the  dimensionless  time  t = vt/X  for 
various  values  of  H . Several  observations  can  be  made;  as  the  value  of  the  Hurst 
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exponent  increases  the  second  spatial  moment  increases  with  dimensionless  time,  a result 
that  is  consistent  with  the  findings  of  Hassan  et  al. 5.  For  each  value  of  the  Hurst  exponent 
the  second  spatial  moment  increases  non-linearly  (non-Fickian)  for  all  simulation  times. 
This  result  can  be  explained  by  the  fact  that  as  the  solute  migrates  in  the  porous  medium, 
it  continuously  encounters  spatially  evolving  scales  of  heterogeneity  (fractal),  hence  the 
dispersion  remains  non-Fickian  for  all  length  scales  below  the  scale  of  heterogeneity  of 
the  modelled  domain  (upper  fractal  cut-off  scale  R ). 


6 Summary 

A two-dimensional  numerical  fractal-based  Monte  Carlo  simulation  model  is  developed  to 
study  flow  and  transport  in  heterogeneous  media  of  spatially  evolving  heterogeneity. 
Realisations  of  the  log  fluctuating  conductivity  field  are  generated  from  a fractal  or  non- 
stationary distribution.  The  long  range  correlation  of  the  hydraulic  conductivity  field  is 
evident  in  the  numerically  obtained  semivariograms.  The  flow  simulations  are  carried  out 
by  solving  a finite  difference  equation  over  a two-dimensional  domain  with  uniform  mean 
flow  in  one  direction.  The  velocity  field  obtained  from  the  conductivity  realisations 
exhibits  long  range  correlations.  A random  walk  particle  tracking  technique  is  used  to 
solve  the  solute  transport  problem  for  non-reactive  solutes.  The  transport  simulations 
indicate  that  higher  Hurst  exponent  H (more  persistence)  of  log  conductivity  results  in 
more  longitudinal  spreading  of  the  contaminant  plume.  Results  also  show  that  irrespective 
of  the  Hurst  exponent  the  transport  process  is  non-Fickian  at  all  times  for  a dispersion 
process  below  the  upper  fractal  cut-off  scale  R . 


References 

1.  Gelhar,  L.  W.  (1993).  Stochastic  Subsurface  Hydrology.  Prentice-Hall,  Englewood 
Cliffs,  New  Jersey. 

2.  Rubin,  Y.  (1990).  Stochastic  analysis  of  macrodispersion  in  heterogeneous  porous 
media.  Water  Resources  Research,  26(1),  133-141. 

3.  Smith,  L.  and  Freeze,  R.  A.  (1979).  Stochastic  analysis  of  steady  state  groundwater 
flow  in  a bounded  domain.  2.  Two-dimensional  simulations.  Water  Resources 
Research,  15(6),  1543-1559. 

4.  Smith,  L.,  and  Schwartz,  F.  W.  (1980).  Mass  transport  1.  A stochastic  analysis  of 
macroscopic  dispersion.  Water  Resources  Research,  16(2),  303-313. 

5.  Hassan,  A.  E.,  Cushman,  J.  H.,  and  Delleur,  J.  W.  (1997).  Monte  Carlo  studies  of 
flow  and  transport  in  fractal  conductivity  fields:  Comparison  with  stochastic 
perturbation  theory,  Water  Resources  Research,  33(11),  2519-2534. 

6.  Beilin,  A.,  Pannone,  M.,  Fiori,  A.,  and  Rinaldo,  A.  (1996).  On  transport  in  porous 
formations  characterised  by  heterogeneity  of  evolving  scales.  Water  Resources 
Research,  32(12),  3485-3496. 

7.  Beilin,  A.,  Salandin,  P.,  and  Rinaldo,  A.  (1992).  Simulation  of  dispersion  in 
heterogeneous  porous  formations:  Statistics,  first-order  theories,  convergence  of 
computations.  Water  Resources  Research,  28(9),  221 1-2227. 


165 


8.  Dagan,  G.  (1990).  Transport  in  heterogeneous  porous  formations:  Spatial  moments, 
ergodicity,  and  effective  dispersion.  Water  Resources  Research,  26(6),  1281-1290. 

9.  Gelhar,  L.  W.,  and  Axness,  C.  W.  (1983).  Three-dimensional  stochastic  analysis  of 
macrodispersion  in  aquifers.  Water  Resources  Research,  19(1),  161-180. 

10.  Sudicky,  E.  A.  (1986).  A natural  gradient  experiment  on  solute  transport  in  a sand 
aquifer:  Spatial  variability  of  hydraulic  conductivity  and  its  role  in  the  dispersion 
process.  Water  Resources  Research,  22(13),  2069-2082. 

11.  Zhan,  H.,  and  Wheatcraft,  S.  W.  (1996).  Macrodispersivity  tensor  for  non-reactive 
solute  transport  in  isotropic  and  anisotropic  fractal  porous  media:  Analytical 
solutions.  Water  Resources  Research,  32(12),  3461-3474. 

12.  Burrough,  P.  A.  (1983).  Multiscale  sources  of  spatial  variation  in  soil.  I.  The 
application  of  fractal  concepts  to  nested  levels  of  soil  variation.  Journal  of  Soil 
Science,  34,  577-597 . 

13.  Molz,  F.  J.,  and  Boman,  G.  K.  (1995).  Further  evidence  of  fractal  structure  in 
hydraulic  conductivity  distributions.  Geophysical  Research  Letters,  22(18),  2545- 
2548. 

14.  Kemblowski,  M.  W.,  and  Wen,  J.  C.  (1993).  Contaminant  spreading  in  stratified  soils 
with  fractal  permeability  distribution.  Water  Resources  Research,  29(2),  419-425. 

15.  Neuman,  S.  P.  (1990).  Universal  scaling  of  hydraulic  conductivities  and  dispersivities 
in  geologic  media.  Water  Resources  Research,  26(8)  1749-1758. 

16.  Mandelbot,  B.  B.  (1983).  The  Fractal  Geometry  of  Nature.  W.  FI.  Freeman,  San 
Francisco. 

17.  Romeu,  R.K.,  and  Noetinger,  B.  (1995).  Calculation  of  intemodal  transmisivities  in 
finite  difference  models  of  flow  in  heterogeneous  porous  media.  Water  Resources 
Research,  31(4),  943-959. 

18.  Chin,  D.  A.,  and  Wang,  T.  (1992).  An  investigation  of  the  validity  of  first-order 
stochastic  dispersion  theories  in  isotropic  porous  media.  Water  Resources  Research, 
28(6),  1531-1542. 

19.  Mantoglou,  A.,  and  Wilson,  J.  L.  (1982).  The  turning  bands  method  for  simulation  of 
random  fields  using  line  generation  by  a spectral  method.  Water  Resources  Research, 
18(5),  1379-1394. 

20.  Yin,  Z.  M.  (1996).  New  methods  for  simulation  of  fractional  Brownian  motion. 
Journal  of  Computational  Physics,  127,  66-72. 

21.  Addison  P.S.  and  Ndumu  A.S.  (1999),  ‘Engineering  Applications  of  Fractional 
Brownian  Motion:  Self-Similar  and  Self-Affine  Random  Processes’,  Fractals,  7(2), 
151-157. 

22.  Addison  P.S.  (1997),  Fractals  and  Chaos:  An  Illustrated  Course,  Institute  of  Physics 
Publishing,  Bristol. 

23.  Tompson,  A.  F.  B.,  and  Gelhar,  L.  W.  (1990).  Numerical  simulation  of  solute 
transport  in  three-dimensional  randomly  heterogeneous  porous  media.  Water 
Resources  Research,  26(10),  2541-2562. 

24.  Voss,  R.  F.  (1989).  Random  fractals:  Self-affinity  in  noise,  music,  mountains,  and 
clouds.  Physica  D,  38,  362-371. 


